/* slfcmbpf.cpp by K.Tsuru */
// function ID 4107 DRADIX since ver 2.182
// first version made on 24 Apr 2010
// remade on 10 Oct 2016 ver 2.30
/****************************************************************
The combination number(binomial function )
           n!
 nCr = ---------
        r!(n-r)!

*****************************************************************/
#ifndef SN_H
#include "sn.h"
#endif
// defined as "const long nMax_combD = 48L;" in "tools.h"
// "const ulong nMax_combL = 2000UL;" in "snmath.h"
SLong comb(ulong n, ulong k){  // combination number(binomial) nCk
    if( n < k ) SNManager::SetError(SNManager::DOMAIN_ERR, "comb()", 4102);
    if(n <= nMax_combD) return (SLong)combD(n, k); // <= 48 uses double vesion
    if(n < nMax_combL) return combL(n,k); // =2000
    return combPF(n,k);
}

/****************************************
PF(Prime factor) method
nCr is evaluated by prime factorization of n!, reduction and product
refer to GMP's factorial algorithm using prime factorization

It makes the prime table between 2 and n
using the sieve of Eratosthenes.
*****************************************/
// for sorting
static int sortByPower(const void *e1, const void *e2) {
    return (int)((const primeFactor *)e2)->power - (int)((const primeFactor *)e1)->power;
}
SLong combPF(const ulong n, const ulong r) {
    if( (r == 0) || (n == r) ) return 1.0;
    if( n < r ) SNManager::SetError(SNManager::DOMAIN_ERR,"combPF()", 4107);

    SNBlock<primeFactor> npf, rpf, n_rpf;
    int in   = FactorialIntoPrimeFactor(n, npf);        // number of prime factor n!
    int ir   = FactorialIntoPrimeFactor(r, rpf);        // r!
    int in_r = FactorialIntoPrimeFactor(n - r, n_rpf);// (n-r)!

    // reduction
    for(int k = 0; k < ir; k++) npf[k].power -= rpf(k).power;
    for(int k = 0; k < in_r; k++) npf[k].power -= n_rpf(k).power;
    int non0pwr = 0;    // number of non zero power elements
    for(int k = 0; k < in; k++) {
        if(npf(k).power) non0pwr++;
    }
    // pick out nonzero power elements
    primeFactor* wpf = new primeFactor[non0pwr+1];
    int el = 0;
    for(int k = 0; k < in; k++) {
        if(npf[k].power) {  wpf[el] = npf(k);   el++;   }
    }
    wpf[el].prime = wpf[el].power = 0;  // index of last element

    qsort(wpf, non0pwr, sizeof(primeFactor), sortByPower);

    SNBlock<primeFactor> pfBk(el+1);

    for(int k = 0; k <= el; k++) pfBk[k] = wpf[k];

    SLong p = PrimeFactorProduct(pfBk, el);
    delete[] wpf;
    return p;
}